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ABSTRACT 


In this study, a soil moisture extraction algorithm was developed 
using a statistical parameter inversion method. Data sets from two air- 
craft experiments were utilized for the test. Both data sets contain 
multi frequency microwave radiometric data as well as surface temperature 
and soil moisture information. 

Results show that by using L, C and X band radiometer data, the 
surface and near surface (<5 cm) soil moisture content can be extracted 
with accuracy of approximately 5-6% for bare fields and fields with grass 
cover. It is demonstrated that this technique is of potential use for 
handling large amounts of remote sensing data from space. 
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1 . INTRODUCTION 


The concept of soil moisture monitoring utilizing passive microwave 
sensors has been investigated both theoretically and experimentally. 
Although both the dielectric properties and radiative transfer processes 
are relatively well understood, the problem still remains complex as 
effects due to various roughness scales, vegetation cover, and other 
mixed terrain features have not been fully modeled and implemented. The 
objective of this study is to take a different approach; instead of 
investigating radiometric responses, feasibility of using a statistical 
inversion scheme for soil moisture retrieval from existing data sets is 
explored. 

The basic tool adopted for retrieving soil moisture information 
from radiometric data is the Statistical Parameter Inversion Method. In 
essence, the method seeks out significant statistical correlations 
between the measured data and their associated parameters. It is a 
general regression technique which minimizes, in the statistical sense, 
the mean square error between the estimated and observed values of the 
parameters of interest. The model is described in Section 2 and the 
computer codes are included in Appendix A. 

Data sets from two aircraft experiments were utilized. Both experi- 
ments were carried out with multi-frequency microwave sensors together 
with detail ground truth measurements. One set of data was collected 
over agriculture fields near Phoenix, Arizona during March 1975. The 
statistical retrieval was carried out for data from bare fields. Radio- 
metric data of 21 cm, 2.8 cm and 1.67 cm, as well as surface temperature 
information, are used to retrieve soil moisture contents in the 0-2 cm 
and 2-5 cm layers. Another set of data was similarly collected over two 
semi-arid areas - Chickasha, Oklahoma and Riesel, Texas - during May 1978. 
The fields were mostly grass-covered. Radiometric data of 21 cm and 5 cm 
as well as surface temperature measurements are used to retrieve volumetric 
soil moisture contents in the 0-2.5 cm and 2.5-5 cm layers. Section 3 
summarizes the data sets and Appendix B illustrates the radiometric 
responses as a function of soil moisture content for the Chickasha/Riesel 
experiment. Further detail of the data can be found in separate reports 
by Burke (1980) and Jackson et al. (1980). 



The results are summarized in Section 4. Over bare fields, the 
soil moisture content can be retrieved with an rms error of approximately 
5% using L and X band measurements as demonstrated from data of the 
Phoenix experiment. Over fields with grass cover, the rms error for 
similar retrieval is approximately 6% using L and C band measurements as 
shown from data of the Chickasha/Riesel experiment. 

In conclusion, the statistical retrieval technique can be valuable 
for handling large amount of data provided the correlation between data 
and parameter is well established either from well controlled experiments 
or theoretical simulation. The technique was demonstrated with two 
aircraft experiment measurements. There are also large amounts of 
multifrequency satellite data available. The next step should be the 
utilization of this type of data for feasibility of soil moisture 
retrieval over large areas. 



2. DESCRIPTION OF THEORY AND MODEL 


In this study, the basic tool adopted for retrieving soil moisture 
information from radiometric data is the Statistical Parameter Inversion 
Method. This technique was originally developed by Westwater and Strand 
(1965) and modified by Waters and Staelin (1968) and Gaut (1967) . It has 
been used extensively in many remote sensing studies and applications. 
Recently, it was employed for development of operational retrievals of 
geophysical parameters ranging from ocean surface wind speed to cloud and 
rain conditions for a seven-channel microwave radiometer (Special Sensor 
Microwave/ Imager (SSM/I)). The SSM/I system will be flown on a future 
mission of the Defense Meteorological Satellite Program. 

In essence, the method seeks out significant statistical correla- 
tions between the measured data and their associated parameters. It is 
a general regression technique which minimizes, in the statistical sense, 
the mean square error between the estimated and observed values of the 
parameters of interest. To do this, it starts from an ensemble of 
simulated cases. These cases are represented by radiometric data 
associated with the physical parameters. The resulting correlations are 
contained in the inversion matrix, which is a set' of coefficients. 

The underlying assumption for this technique is that there exist 
some linear combination of data elements which will provide a useful 
estimate of the desired parameters. That is to say: 


^ U 

p. = y D..d. 
^ i=i J 


where p^ is the estimate of the parameter p^, 
temperatures and D. . is the inversion matrix, 
vectorial ly as: 


dj's are the brightness 
This can be written 


£ = D • d 

The determining conditions on the inversion matrices results from 
minimizing the square of the deviation of the estimated from the actual 
values of the parameters over the entire ensemble of simulated scenes. 
It can be shown as: 

D = C(p,d) • C"^ (d,d) 



where C(p,^) and are the correlation matrices defined as 

C- . (x,y) = <x. ,y.> 

As can be seen, the Statistical Parameter Inversion Method relies on 
^ priori information defining the correlation between the predicted data 
vector and the predicted geophysical parameter vector (£) • 

Another feature can be used in the inversion process to "tune" the 
inversion matrix to the average value of each parameter. To do this, 
the data vector is transformed like: 

I (d) = (1, d^ - dj, - d^, .... d^ - d^) 

The effect of this transformation is to add a column to the inversion 
matrix which contains the ensemble averages of the parameters. Though 
the resulting inversion is still linear in the T 's, D is now focused on 

D = 

correlated deviations from the expected "mean" values of the parameters. 
For example, in retrieving soil moisture, the expected mean value can be 
that usually occurring or expected for the season. However, data for 
this study were acquired from a number of experiments which included a 
broad range of soil moisture conditions. Therefore, the "tuning" 
process to a mean or most probable case was not used. 

Problems arise which can degrade the inversion results when there 
are non-linearities in the ^ to £ relationships. Various modifications 
can be developed to correct for this problem (e.g. Staelin, 1980; Wilheit 
and Chang, 1980). In the remote sensing of soil moisture in the passive 
microwave region, the brightness temperature to soil moisture relation- 
ship is approximately linear; therefore, treatments for non-linearities 
were not considered for this study either. 

Appendix A lists the FORTRAN codes of the computer program developed 
for this study. It consists of five subroutines with their main functions 
listed as follows: 

1) MAIN-Program SPIM: the Statistical Parameter Inversion 

Method is used to obtain the inversion matrix, upon 
option, results are then evaluated by computing the 
standard deviation between estimated (computed) and 
actual parameter values; 



2) INDATA: to set up data and parameters for inversion 

matrix calculation; 

3) CORMET: to calculate the correlation matrices from 

the known parameter and data sets; 

4) MINV: an adopted IBM application program to perform 

matrix inversion using the standard Gauss-Jordan 
method; and 

5) GMPRD: to multiply two matrices to form a resultant 

general matrix. 



3. DESCRIPTION OF DATA 


In this study, data from two aircraft experiments were utilized. 

One was carried out for agriculture fields near Phoenix, Arizona during 
March 1975. Another set of data was collected in two semi-arid areas - 
Chicasha, Oklahoma and Riesel, Texas during May 1978. Data from the 
Phoenix experiment was described and analyzed in a separate report by 
Burke (1980). Data from the Chickasha and Riesel experiment was 
described in Jackson et al. (1980). The following is a summary of the 
two data sets. 

3.1 Data from the Phoenix Experiment 

During March 1975, an aircraft mission consisting of four flights 
over the Phoenix, Arizona test site was conducted for the purpose of 
studying the use of microwave radiometers for the remote sensing of soil 
moisture . The investigators involved in this mission came from NASA, 
the Agricultural Research Service of USDA, the University of Arkansas, 
and Texas ASM University, This mission consisted of predawn and midday 
flights on 18 and 22 March 1975. There were radiometers operating at 
wavelengths of 0.8, 1.3, 1.67, 2.8 and 21 cm. The 2.8 cm radiometer is 
a dual -polarized conically scanning radiometer operating at a fixed look 
angle of 50°. The other radiometers were non-scanning but could have 
their nadir look angles varied. In addition to the microwave instruments 
the scientific package included a thermal infrared radiometer for measur- 
ing surface temperature, Three passes were taken over each field; one at 
a nadir angle of 0 = 0° and two at a nadir angle of 0 = 40° alternating 
the polarization sensitivity of the 21 cm antenna. 

Ground measurements were made in 46 fields. Twenty-eight were 
without vegetative cover and 18 had vegetative covers of either alfalfa 
or wheat. The fields, which have an area of 16 hectares (400 x 400 m) , 
were arranged in pairs to provide a uniform target 800 meters wide. The 
soil moisture sampling procedures for this mission included measurements 
of the moisture content and temperature in each of the following soil 
layers: 0-1 cm, 1-2 cm, 2-5 cm, 5-9 cm, and 9-15 cm. The results 
included data from a variety of moisture conditions due to the irrigation 
and drying cycles of the fields. 



The data analysis was carried out for 21 cm, 2.8 cm and 1.67 cm 
measurements. Analysis of data from 1.3 cm and 0.8 cm were not carried 
out as they were relatively insensitive to soil moisture content. 

Responses to both bare and vegetated fields for nadir and polarized 
conditions were investigated in detail in Burke (1980) . 

In this study, soil moisture contents in the 0-2 cm and 2-5 cm 
layers are taken to be the two parameters to be retrieved from statistical 
inversion; they represent the surface and near surface moisture conditions 
Only bare fields are treated. The average physical temperature of the 
surface is used to normalize the brightness temperature such that morning 
and afternoon measurements can be utilized simultaneously. 

3.2 Data from the Chicasha and Riesel Experiment 

Cooperative investigations were conducted during May 1978 by the 
National Aeronautics and Space Administration (NASA) and the U.S. Depart- 
ment of Agriculture (USDA) as part of a project to evaluate remote 
sensing in hydrologic studies with primary emphasis on measurements. 

Ground observations and aircraft remote sensing experiments were con- 
ducted in two semi-arid areas - Chicasha, Oklahoma and Riesel, Texas. 
Participants in the study were the NASA Goddard Space Flight Center and 
the USDA-SEA-AR Hydrology Laboratory, Southern Plains Watershed and 
Water Quality Laboratory and the Grassland, Soil, and Water Laboratory. 

Three successful flights were made on May 1, 12 and 30, over the two 
test sites. The study sites and watersheds in Oklahoma were in the 
Washita River Experimental Watershed, Chickasha, area. Most watersheds 
had a dense grass cover; some were bare or wheat fields. The study site 
located in the central part of Texas near Riesel was also on an experi- 
mental watershed area. ^ The land cover varied from almost bare soil to a 
very dense vegetative cover. Some fields experienced considerable 
vegetation growth between the first flight and the last flight. 

Sensor data included those from photography, thermal infrared radio- 
meter, passive microwave sensor system of L-band (21 cm), C-band (5 cm), 
X-band (2.8 cm), Ku-band (1.67 cm), K-band (1.35 cm) and Ka-band (.81 cm) 
with look angles at 0° and/or 40°, and scatterometer operating at L, C 
and P bands. For fields from both study sites, all soil moisture 



samples were weighed, oven-dried and weighed again to determine their 
gravimetric soil moisture. Bulk density samples from each field were 
then used to determine the volumetric soil moisture. Soil moisture 
information in the 0-2.5 cm, 2.5-5 cm, 5-10 cm and 10-15 cm was obtained. 
(Jackson et al . , 1980). 

In this study, volumetric soil moisture contents in the 0-2.5 cm 
and 2.5-5 cm are taken as the two parameters to be retrieved. Surface 
temperature obtained from the infrared sensor is used to normalize 
different conditions from various days and sites. Radiometric data from 
L, C, X and Ku bands are utilized. Since there was no previous study on 
analysis of the data, brightness temperatures at these wavelengths are 
plotted as a function of the 0-2.5 cm layer volumetric moisture content 
(Appendix A) for a general feeling of the radiometric response. 



4. RESULTS AND CONCLUSIONS 


Figui-e 4-1 demonstrates an ideal approach of parameter (e.g. soil 
moisture) retrieval from remote sensing data. An ensemble of the physical 
system is set up first. This system contains a broad and representative 
range of model cases, each with a different combination of soil moisture 
condition, surface roughness, vegetation coverage, temperature, atmo- 
spheric contamination and other factors that can affect the radiometric 
signature. Applying the radiative transfer and other model computations, 
an ensemble of simulated radiometric data set is obtained. Correlation 
matrices of soil moisture and radiometric data as well as the inversion 
matrix are calculated as outlined in Section 2. Measured radiometric data 
set is then applied to the inversion matrix to obtain the estimated 
parameter values. The estimated values can then be compared with the 
ground truth data for error analysis. 

In the present study, the approach is somewhat different due to the 
following; 

1) Data from well controlled experiments are utilized. 

The data amount is relatively small; variability in 
soil moisture is broad, but much less significant in 
other conditions within the same experiment; for 
example, soil types and physical temperatures. It is 
thus not applicable to set up a statistical ensemble 
with model computation of radiometric prediction for 

the inversion matrix. 

2) The main purpose of the study is to test the retriev- 
ability of soil moisture condition from radiometric 
data rather than to test the accuracies of theoreti- 
cal modeling. Theoretical studies of soil moisture 
response have been carried out in great extent. In 
Burke (1980) report, data from the Phoenix experiment 
showed good agreement with model computation. 

As a result, measured data are used both for the derivation of the inver- 

2 

sion matrix and also error analysis. The rms error defined by [Z(p-p*) /N] 
is used where p and p* are the parameter and the estimated value and N is 
the number of cases. 




Figure 4-1 Flow diagram for parameter retrieval using Statistical 
Parameter Inversion Method 









4.1 Results of Soil Moisture Retrieval Utilizing Data 

from the Phoenix Experiment 

As described in Section 3.1, soil moisture contents in the 0-2 cm 
and 2-5 cm layers are taken as the two parameters to be retrieved. 

Tests are carried out using (1] five channels of data from 21 cm (0°, 
vertical and horizontal) and 2.8 cm (vertical and horizontal) and ( 2 ) 
eight channels with additional data from 1.67 cm (0°, vertical and hori- 
zontal). Since emissions at 2.8 cm and 1.67 cm from vegetated fields are 
relatively insensitive to soil moisture contents, retrievals are carried 
out for bare fields only. 

Question arises as to whether the data sets for deriving the 
inversion matrix and for accuracy test can be the same or they should 
be independent of each other. The morning measurements from two separate 
days are used for demonstration. Table 1 shows the difference between 
using data separately and combined. Data from five channels are 
utilized in this case. As can be seen, using the same data sets for 
derivation and test causes the rms error of retrieval to reduce by 
1 - 1 ^ 5 % in soil moisture content. Bearing this in mind, the rest of the 
retrievals are carried out using identical sets of data for both deriva- 
tion and accuracy test so that maximum number of cases and range of data 
are used. 


TABLE 1 

TEST OF RMS INVERSION ERRORS OF SOIL MOISTURE CONTENT FOR 
VARIOUS WAYS OF USING DATA SETS FOR SET UP AND EVALUATION 


Data for 
Inversion Matrix 
Set up 

Data for 
Accuracy 
Test 

RMS Inversion Error 
(% Soil Moisture) 

0-2 cm 

2-5 cm 

Mar 18 AM (14)* 

Mar 22 AM (13) 

5.6 

5.1 

Mar 22 AM (13) 

Mar 18 AM (14) 

Mar 18 AM 
Mar 22 AM 

Mar 18 AM ^^ 7 ) 
Mar 22 AM 


4.0 


*number in parenthesis refers to the number of cases 


Table 2 demonstrates the utilization of both morning and afternoon 
measurements with five and seven channels of data. Two approaches are 
tested. In the first case, the brightness temperatures are directly 
used. In another case, the brightness temperatures are first normalized 
















to the average surface temperatures; 278°K, 300°K, 280‘’K and 296°K are 
used for March 18 AM, PM and March 22 AM, PM measurements, respectively. 
Fifty-five cases are available. For five channels, the rms errors for 
0-2 cm and 2-5 cm moisture contents are 5.1% and 4.6%, respectively; a 
degradation as compared to the 4.3% and 4.0% for morning measurements 
alone due to the differences in moisture distribution in soil between 
mornings and afternoons. The improvements from normalizing the bright- 
ness temperatures to surface temperatures (T^) are minimal, from 5.1% and 
.4.6% to 5.0% and 4.3%. 


TABLE 2 

RETRIEVAL ERROR OF MOISTURE CONTENT BY WEIGHT IN % 
FOR THE PHOENIX EXPERIMENT, 1975 


Radiometric Data* 


RMS Inversion Error 

0-2 cm 

2-5 cm 

Tb 

Tb/Ts 

Tb 

Tb/Ts 

21, 2.8 cm (5) 

Mar 18, 22 
AM, PM 

■ 

5.0 

B 

B 

21, 2.8, 1.67 cm 
C8) 

Mar 18, 22 
AM, PM 

B 

m 


B 


*number in parenthesis refers to the number of channels 
**number in parenthesis refers to the number of cases 


Results also show that by utilizing the additional 1.67 cm data, the 
0-2 cm error is reduced by about 1^% and had little impact on the 2-5 cm 
results due to its sensitivity to surface moisture content only. The 
2-5 cm results have been consistently somewhat better than those of the 
0-2 cm. This is due to the less variability in the subsurface moisture 
conditions between different fields and measurements. Although radio- 
metric data are not more sensitive to subsurface moisture conditions, 
the statistical inversion error is nevertheless smaller. 

4.2 Results of Soil Moisture Retrieval Utilizing Data from 
the Chickasha/Riesel Experiment 

Similar procedure is adopted for the Chickasha/Riesel data set. 
Volumetric soil moisture contents in the 0-2.5 cm and 2.5-5 cm layers 
are the two parameters to be retrieved. From the data presented in 




















Appendix B, it is noted that at 1.67 and 2.8 cm the amount of missing 
data is high and the sensitivity of the available data to surface moisture 
is also relatively poor. Data from the five channels of 21 cm (0° and 
horizontal) and 6 cm (0®, vertical and horizontal) (0 = 40°) are thus util- 
ized. There are cases with missing data in some of the five channels. 

Since both the derivation of the inversion matrix and also the parameter 
estimation do not allow missing data, adjustments have to be made. It is 
decided that if data from more than two out of the five channels are 
missing, the case is discarded; otherwise, the missing brightness temp- 
erature is replaced by the average measured value of that area (Chickasha 
or Riesel) on that day. As a result, 41 cases are available. As shown 
in Table 3, the rms inversion errors of the 0-2.5 cm, 2.5-5 cm layers 
are 7.5% and 6.3%, respectively. 


TABLE 3 

RETRIEVAL ERROR OF VOLUMETRIC MOISTURE CONTENT IN % 
FOR THE CHICKASHA/RIESEL EXPERIMENT, 1978 


Radiometric 

Data* 

Data Set** 

RMS Inversion Error (%) 

0-2.5 cm 

2.5-5 cm 

21, 6 cm 
(5) Tg's 

May 1, 12, 30 
(41) 

7.5 

6.3 

21, 6 cm 
(5) Tg's 

May 1, 12, 30 
(39) 

6.1 

5.4 

21, 6 cm 

(5)(Tg/T<,)'s 

May 1, 12, 30 
(39) 

6.0 

5.5 


*number in parenthesis refers to the number of channels 
**number in parenthesis refers to the number of cases 


A closer look of the plotted data show that the two fields with 
highest moisture contents (53,7% and 50.9% for RG 83, RG 88, respectively 
on May 30) are totally uncorrelated to the rest of the data set. The 
reason is not known. These two cases are removed and the new retrieval 
showed rms inversion errors of 6,1% and 5.4% for moisture contents in 
the 0-2.5 cm and 2.5-5 cm layers, respectively (Table 3). 

The same 39 cases are tested with the brightness temperatures norma- 
lized to their surface temperature (T^) measurements. The results are 
also shown in Table 3 with inversion errors of 6.0% and 5.5%, virtually 
the same as the previous case. 






4.3 Further Discussion 


In Sections 4.1 and 4.2, the results of retrieval errors from two 

aircraft experiments were presented. One major concern is that by using 

the same data set for both deriving the algorithm and testing it, the 

uncertainty levels are just those associated with the noise in the data 

set and do not indicate the uncertainty introduced by the extraction 

algorithm. Therefore, one more set of tests using the Phoenix data was 

carried out in order to delineate the difference between uncertainties 

from the extraction algorithm and the noise in the data. Furthermore, 

2 % 

the relative rms error, defined as [ICp-p*/p) /N] > was also computed. 

The results are shown in Table 4 using normalized brightness temperatures 
from all eight channels. 


TABLE 4 

RMS ERROR, BOTH ABSOLUTE AND RELATIVE, 

OF SOIL MOISTURE RETRIEVAL, OF THE PHOENIX EXPERIMENT, 1975 


Data for 
Inversion 
Matrix Setup 

Data for 
Accuracy 
Test 

RMS Inversion Error 

WEmssasmm 

Relative 

0-2 cm 

2-5 cm 

0-2 cm 

2-5 cm 

Mar 18, 22 (55)* 
AM, PM 

Mar 18, 22 (55) 
AM, PM 

m 

B 

.41 

.27 

Mar 18 
AM, PM (28) 

Mar 22 
AM, PM (27) 

4.7 

m 

.50 

.35 


*number in parenthesis refers to the number of cases 


As can be seen again, the main source of inversion error is from 
the noise within the data set rather than the extraction algorithm. It 
also demonstrated that data sets from the two different days (March 18 
and March 22) are compatable with each other. The relative rms errors 
are generally within 50% (.5) and better for the subsurface retrievals 
('\'30% or .3). Figure 4-2 is a scatter plot of calculated soil moisture 
values versus those measured for the second test in Table 4. 























CALCULATED SOIL MOISTURE CONTENT (%) 
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Figure 4-2 




4.4 Conclusions 


In this study, the retrievability of soil moisture content from 
remotely sensed radiometric data is demonstrated. Over bare fields, the 
moisture content can be retrieved with an rms error of approximately 5% 
using L (21 cm) and X (2.8 cm) band measurements as shown from data of 
the Phoenix experiment. Over fields with vegetation cover, the rms error 
for similar moisture retrieval is approximately 6% using L (21 cm) and 
C (6 cm) band measurements as shown from data from the Chickasha/Riesel 
experiment. The results are promising as actual field data are used and 
only unbiased mathetmatical inversion is applied. 

One unexpected result is the lack of improvement when surface 
temperature measurements are incorporated to normalize the brightness 
temperatures. This is probably due to the fact that at the wavelengths 
considered, the physical temperature factor contributed to the signature 
is an integrated surface-subsurface component, especially at X and C band 
wavelengths, rather than the surface temperature alone. As all measure- 
ments were taken with relatively uniform subsurface temperature, the con- 
tribution of the physical temperature may well be within noise level as 
of other factors (roughness, soil type, etc.). 

The handling of missing data from some of the channels is important 
as it is an inevitable problem of inversion from measured data. In this 
study, average value of data of similar time and space domains is used. 
This can be easily adopted for operational use. 

Another restriction of the study is the limited amount of useful 
data. In the Phoenix data set, 55 cases are available; and for the 
Chickasha/Riesel data set, 39. Ideally, there should be independent 
data sets for the derivation of inversion and estimation of parameters. 

An attempt was initiated to correlate the two data sets. It was quickly 
recognized to be not feasible due to the differences between the two 
data sets in the wavelengths used, the nature of background (bare versus 
grassland), and ground truth measurements (layer thicknesses, moisture 
contents by weight versus volume). 

In conclusion, the technique presented and tested in this study 
proves to be of potential use for large amount of data retrieval. There 
have been data of this nature from spacebome sensors; for example, SMMR 



(Scanning Multichannel Microwave Radiometer) from both Seasat and Nimbus-7 
satellites. After initial analysis of the data, such inversion technique 
for soil moisture retrievals, as well as for other parameters, should be 
implemented. 
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APPENDIX A 


FORTRAN CODES OF THE 
STATISTICAL PARAMETER INVERSION METHOD 



oonooooooonooonooooooooooooonnoooooooooooooooo 


■CGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
PROGRAM SPIM 
PURPOSE: 

STATISTICAL PARAMETER INVERSION METHOD IS USED TO 
OBTAIN THE INVERSION MATRIX AND PARAMETER ESTIMATIONS 

DISCRIPTIONS OF PARAMETERS: 

NCASE - NUMBER OF KNOWN CASES (MAX r 500, INPUT) 

NP - NUMBER OF PARAMETERS (MAX = 10, INPUT) 

ND - NUMBER OF DATA IN A CASE (MAX = 10, INPUT) 

DFLAG - .TRUE. FOR PARAMETER ESTIMATION, .FALSE. FOR MATRIX 
ESTIMATION ONLY 

NINP - NUMBER OF CASES FOR PARAMETER ESTIMATION 
EFLAG - .TRUE. FOR MATRIX EVALUATION, .FALSE. OTHERWISE 
PARA - ARRAY OF PARAMETER 

DATA - ARRAY OF DATA 

CORPO - CORRELATION ARRAY OF PARAMETER AND DATA 
CORDD - CORRELATION ARRAY OF DATA AMD DATA 
DMT - INVERSION MATRIX 
DET - DETERMINANT OF MATRIX INVERSION 
LWOK - TEMPORARY WORKING ARRAY FOR MINV 

MWOK - TEMPORARY WORKING ARRAY FOR MINV 

PEVA - INPUT parameter ARRAY FOR EVALUATION 
OIFF - DIFFERENCE BETWEEN ESTIMATED AND INPUT PARAMETER 
VALUES 

SDV - STANDARD DEVIATION OF EVALUATION 

SUBROUTINE AND FUNCTION SUBPROGRAMS: 

INOATA - ROUTINE TO INPUT NCASES OF DATA AND PARAMETERS 
COREMT - ROUTINE TO CALCULATE CORRELATION MATRICES 
MINV - ROUTINE TO CALCULATE MATRIX INVERSION 
GMPRD - MATRICES MULTIPLICATION 

REMARKS: 

INPUT NCASE, NP, NO, DFLAG, NINP, EFLAG IN F0RMAT(/ai5,L2, I5,L2) 
METHOD 

STATISCAL PARAMETER INVERSION TECHNIQUE IS USED TO GET 
INVERSION MATRIX, RESULTS ARE EVALUATED BY‘ COMPUTING 
STANDARD DEVIATION. 

DESIGNER: H.K, BURKE AND J-H HO 

PROGRAMMER: J-H HO 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
logical DFLAG, EFLAG 
INTEGER NP, NO, NCASE, NINP 
DIMENSION PARA(5000),DAT(5000) 

DIMENSION CORPO (100), CORDD (100), DMT (100) 

DIMENSION LWOK(10),MWOK(10) 

DIMENSION DIFF(10),PEVA(10),SDV(10) 

DATA SDV/10*0./ 

C 

C 

READ(5,5000) 

5000 FOKMAT(IOX) 

READ (5, 5001) NC ASE , NP, NO, DFLAG, NINP , EFLAG 

5001 F0RMAT(3I5,L2,I5,L2) 

WRITE (t), 6001) NCASE, NP,ND, DFLAG, NINP, EFLAG 
6001 FORMAT(1H1,OX,6HINPUT:,/,10X,6HNCASE=,I5,/, 10X,3HNP=,I8, 



6000 

6002 

20 


5002 

ao 


6003 

6004 


50 

6005 

6006 
60 

60 

6007 

100 


• /a0X,3HNO=,I8,/»10X,6HOFLAG=,L5»/#10X,5HNINP=,I6, 

/#10X,6HEFLAG=,L5) 

READ(5#5000) 

call iNDATA(PARA,NP,OAT,NO,NCASE) 

CALL COREMT(PARA/NP,DAT,NO,NCASE,CORPD,COROD) 

CALL MINV(CORDO,ND,D£T#LWOK,MWOK) 

if(det,eq,o) stop 

CALL GMPRD(CORPOfCOROD,DMT»NP,ND,NO) 

IEND=(ND-1)*NP 

WRITE(6,6000) 

F0RMAT(//2X,16HMATRIX ELEMENTS:) 

DO 20 1=1, NP 

IENOsIENO+1 

WRITE (6, 6002) (DMT ( J) , Jsl , lEND, NP) 

F0RMAT(2X, lOEll.4) 

CONTINUE 

IF(,NOT,OFLAG) go to 100 
READ(5,5000) 

IEND=0 

DO 40 1=1,NINP 
ITEMsIEND+1 
IEN0=IEN0+N0 

READ (5, 5002) (DAT ( J) , J=ITEM, lEND) 

F0RMAT(10F8.3) 

CONTINUE 

CALL GMPRD(DMT,DAT,PARA,NP,ND,NINP) 

IF(EFLAG) READ(5,5000) 

lENOsO 

IEN02=0 

DO 60 I=1,NINP 
iTEMsIEND+l 
IEN0*IEND+N0 
ITEM2SIEND2+1 
IEND2=IEND2tNP 

WRITE (6, 6003) (DAT( J) , J=ITEM, lEND) 

FORMAT (//4X, 11 HINPUT DATA : , /2X , 1 OF 1 1 , 4) 

WRITE (6,6004) (PARA( J) , J=ITEM2, IEND2) 

FORMAT (4X,12H0UTPUT PARA : , /2X , 1 OF 1 1 ,4) 
IF(.NOT.EFLAG) GO TO 60 

READ(5,5002) (PEVA ( J) , Jb 1 , NP ) 

DO 50 Kb1?NP ' - ' - 

0IFF(K)=PARA(ITEM2+K-1)-PEVA(K) 

SDV(K)bSDV(K)+OIFF(K)**2 

CONTINUE 

WRITE(6,60D5) (PEVA ( J) , J= 1 , NP ) 

FORMAT (4X,12HKN0WN PARA : , /2X, 1 OFl 1 .4 ) 
WRITE(6,6006) (DIFF( J) , J=1 ,NP) 

FORMAT (2X, 5HDI FF : ,46X,10F6. 2) 

CONTINUE 

IF(,NOT,EFLAG) GO TO 100 
DO 80 Ksl,NP 

SDV(K)=SORT(SDV(K)/NINP) 

CONTINUE 

WRITE(6,6007) (SDV (K) ,Ksl ,NP) 

FORMAT(///,22HEVALUATION SDV VALUES: ,/2X, lOFl 1 ,4) 
CONTINUE 

STOP 

END 
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SUBROUTINE INOATA (PARA , NP, DAT, NO r NCASE) 

ccccdcccctcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE INOATA 

PURPOSE: 

INPUT DATA AND PARAMETER FOR INVERSION MATRIX CALCULATION 

USAGE: CALL INDATA (PAR A , NP, DAT, NO , NCASE ) 

OISCRIPTION OF PARAMETERS: 

PARA - ARRAY OF PARAMETERS (OUTPUT) 

DAT- ARRAY OF DATA (OUTPUT) 

NCASE - number of CASES (INPUT) 

NP - NUMBER OF PARAMETERS FOR A CASE (INPUT) 

NO - NUMBER OF DATA FOR A CASE (INPUT) 

SUBROUTINE AND FUNCTION SUBPROGRAMS: 

NONE 


REMARKS: 

, INPUT A CASE AT A TIME, USING FORMAT ( I 0F8 , 3) FOR DATA FIRST, 

' THEN INPUT PARAMETER BY FORMAT ( 1 OFB , 3) NEXT LINE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DIMENSION PARA(NP, NCASE), OAT (NO, NCASE) 

WRITE(6,600l) 

6001 F0RMAT(1H1,5X, 12HINPUT CASE :/) 

DO 20 Isl, NCASE 

READ(5,500l) ( DAT ( J, I ) , Jsl , ND) 

/ REAO(5,5001) (PARA(J,I),Js1,NP) 

5001 FORMAT(IOF8,3) 

WRITE(6,6002) (DAT ( J, I ) , Jsl , ND) 

WRITE(6,6003) (PARA ( J, I ) , J=1 , NP) 

6002 FORMAT(/2X,AHDAT:,/2X,10Fll.a) 

6003 FORMAT(2X,5HPARA:,/2X,lOFll.a) 

20 CONTINUE 

RETURN 

END 
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- SUBROUTINE COPEMT (PARA , NP, DAT, ND, NCASE , CORPD/ COROD) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE COREMT 

PURPOSE 

TO CALAULATE THE CORRELATION MATRICES FROM 'NCASE' SET OF 
KNOWN PARAMETERS AND DATA 

USAGE 

CALL COREMT(PARA,NP,OAT,ND,NCASE,CORPD,CORDD) 

DISCRIPTION OF PARAMETERS: 

PARA - ARRAY OF PARAMETERS 
NP - NUMBER OF PARAMETERS FOR A CASE 
DAT - ARRAY OF DATA 
ND -NUMBER OF DATA FOR A CASE 
NCASE - NUMBER OF CASES OF DATA 
CORPD - CORRELATION ARRAY OF PARAMETER AND DATA 
CORDO - CORRELATION ARRAY OF DATA AND DATA 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED: NONE 

REMARKS! NONE 

METHOD: 

CORPDd, J)=SUM OF PARA(I,K)*DAT(J,K) OVER K (NCASE) 

CORPDd/ J)s8UM OF OAT (I , K) *DAT ( J, K ) OVER K (NCASE) 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DIMENSION PARA(NP,NCASE),DAT(NDr NCASE) 

DIMENSION CORPD ( NP, NO) , CORDO (NO, NO) 

C INITIALIZATION 

DO 50 J=1,ND 
DO 10 lsl,NP 

CORPDd, J)=0. 

10 CONTINUE 

DO 20 Isl,ND 

CORDOd, J)sO. 

20 CONTINUE 

50 CONTINUE 

C 

DO 200 K=l, NCASE 
DO 150 J=l,NO 

DO 120 Isl,NP 

CORPD(I,J)*CORPD(I, J)+PARA(I,K)*DAT(J,K) 

120 CONTINUE 

DO 140 1=1, ND 

CORDOd, J)sCOROO(I, J)+DAT( I, K)*DAT(J,K) 

140 CONTINUE 

150 CONTINUE 

200 CONTINUE 

C 

DO 300 J=l,ND 
DO 220 1=1, NP 

CORPDd, J)=CORPD( I, J) /NCASE 
220 CONTINUE 

DO 240 1=1, ND 

CORDD(I,J)=CORDD(I,J) /NCASE 
240 CONTINUE 

300 CONTINUE 

RETURN 
END 
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• SUBROUtiNE MINV(ArN,D,L»M) 
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NAME: MATRIX INVERSE UTILITY (MINV) 

PURPOSE: TO PERFORM MATRIX INVERSION USING THE STANDARD 

GAUSS-JORDAN METHOD. THE DETERMINANT IS ALSO CAL- 
CULATED. A DETERMINANT OF ZERO INDICATES THAT THE 
MATRIX IS SINGULAR. 


INTERFACES: CALLING MODULES SHOULD HAVE A, L, M PROPERLY 

DIMENSIONED. 


CALLED MODULES - NONE 


INPUT PARAMETERS - A = INPUT MATRIX 

N = ORDER OF MATRIX A 

OUTPUT PARAMETERS - D s RESULTANT DETERMINANT 

L = WORKING MATRIX 
M = WORKING MATRIX 
A s OUTPUT MATRIX 


COMMON BLOCKS SET - NONE 


COMMON BLOCKS READ - NONE 
DATA FILES - NONE 


COMMENTS: IBM APPLICATION PROGRAM 

SYSTEM/560 SCIENTIFIC SUBROUTINE PACKAGE 
(360A-CM-03X) VERSION III 
PROGRAMMER'S MANUAL 


*********************************************************************** 


DIMENSION A(n,L(l),M(l) 


( 


{ 


( 

I 


c. 




c ********** SEARCH FOR LARGEST ELE^MENT 

D = 1.0 
NK = -N 
DO 195 K=1,N 
NK r NK+N 
L(K) = K 
M(K) = K 
KK s NK+K 
6IGA s A(KK) 

DO 115 J=K,N 
IZ = N*(J-1) 

DO 115 I=K,N 
IJ = IZ+L 

105 IF (ABS(BIGA) - ABS(A<IJ))) 110/115»115 
110 8IGA * A(IJ) 

L(K) = I 
M(K) = J 
115 CONTINUE 

C ********** INTERCHANGE ROWS 


***** 


***** 


V 


J = L(K) 

IF (J-K) 130,130»120 



4 


. .. . 120. KI 5 K-N 

00 125 IP1,N 
KI s KI+N 
HOLD = -A(KI) 
f JI s KI-K+J 

A(KI) P A(JI) 

125 A(JI) p HOLD 
r 

C ********** INTERCHANGE COLUMNS 

^ 130 I = M(K) 

IF (I-K) t«5fia5,135 
135 JP s N*(I-1) 

( DO lao JP1,N 

JK = NK+J 

JI P JP+J 

< HOLD = -A(JK) 

A(JK) P A(JI) 
lAO A(JI) = HOLD 

f 

C ********** DIVIDE- COLUMN BY MINUS PIVOT(BIGA) 

{ 1A5 IF (ABSCBIGA-.OOOU.GT.O,) GO TO 155 

150 D = 0,0 

WRITEC6, 1000) 

f 1000 F0RMAT(///,5X, 16H SINGULAR MATRIX//5X, 12H DETERMINANT, 

+ 17H LESS THAN 0,0001) 

RETURN 

C 155 DO 165 islrN 

IF (I-K) 160,165,160 
160 IK P NK+I 

C A(IK) p A(IK)/(-BIGA) 

165 CONTINUE 

i C ********** REDUCE MATRIX 

DO 180 1=1, N 
I IK = NK+I 

HOLD = A(IK) 

IJ = I-N . 

i DO 180 J=1,N 

IJ = IJ+N 

IF (I-K) 170,180,170 
{ 170 IF (J-K) 175,180,175 

175 KJ = IJ-I+K 

A(IJ) = HOLO*A(KJ)+A(IJ) 
i_ 180 CONTINUE 

C ********** DIVIDE ROW BY PIVOT 

L 

KJ = K-N 
DO 190 J=1,N 

V KJ = KJ+N 

IF (J-K) 185,190,185 
185 A(KJ) = A(KJ)/BI6A 

V 190 CONTINUE 

C ********** PRODUCT OF PIVOTS 

D = D*BIGA 


***** 


***** 


***** 


***** 


***** 


V, 


C ********** 


REPLACE PIVOT BY RECIPROCAL 


***** 



. A(KjO = 1,0/BIGA 
195 CO.MTINUE 

********** FINAL ROW AND COLUMN INTERCHANGE ***** 

K = N 

200 K = (K-n 

IF (K) 100,100,205 
205 I = L(K) 

IF (I-K) 220,220,210 
210 JO = N*(K-1) 

JR = N*(I-1) 

DO 215 J=1,N 
JK = JO+J 
HOLD = A(JK) 

JI = JR+J 
A(JK) = -A(JI) 

215 A(JI) = HOLD 
220 J = M(h) 

IF (J-K) 200,200,225 
225 KI = K-N 

DO 230 1=1, N 
KI = KI+N 
HOLD s A(KI) 

JI = KI-K+J 
A(KI) = -A(JI) 

230 A(JI) = HOLD 
GO TO 200 
100 RETURN 
END 
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. SUBROUTINE GMPRD ( A , B» R , N, M,L) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
SUBROUTINE GMPRD 
PURPOSE 

MULTIPLY TWO MATRICES TO FORM A RESULTANT GENERAL MATRIX 
USAGE 

CALL GMPRD(A,8,R,N,M,L) 

DESCRIPTION OF PARAMETERS 

A - NAME OF FIRST INPUT MATRIX 
B - NAME OF SECOND INPUT MATRIX 
R - NAME OF OUTPUT MATRIX 
N - NUMBER OF ROWS IN A 

M • NUMBER OF COLUMNS IN A AND ROWS IN B 
N - NUMBER OF COLUMNS IN B 


REMARKS 

ALL MATRICES MUST BE STORED AS GENERAL MATRICES 

MATRIX R CANNOT BE N THE SAME LOCATION AS MATRIX A 

MATRIX R CANNOT BE N THE SAME LOCATION AS MATRIX B 

NUMBER OF COLUMNS OF MATRIX A MUST BE EQUAL TO NUMBER OF 

ROWS OF MATRIX B 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

THE M BY L MATRIX 8 IS PREMULTIPLIED BY THE N BY M MATRIX 
A AND THE RESULT IS STORED IN THE N BY L MATRIX R, 

SYSTEMS 

IBM/360 IN SCIENTIFIC SUBROUTINE PACKAGE 
ADOPTED BY J-H HO, 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c , ^ 

c 

DIMENSION A(l),B(U,Ra) 

IRsO 

IKs-M 

DO 10 Kcl,L 

IKsIK+M 

DO 10 Jsl,N 

IRsIR +1 

JI*J»N 

IBsIK 

R(IR)=0 

DO 10 lBl,M 

JIsJUN 

IBbIB+1 

10 R(IR)=R(IR)+A(JI)*B(IB) 

RETURN 

END 



APPENDIX B 


RADIOMETRIC RESPONSES AS A 
FUNCTION OF SURFACE SOIL MOISTURE CONTENT 
(Chickasha/Riesel, 1978) 
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VOLUMETRIC SOIL MOISTURE CONTENT (%) 

Figure B-1 The Ku band (1.67 cm) radiometer data (look angle = 0° and 
40°) as a function of volumetric soil moisture content (%) 
in the 0-2.5 cm layer 
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Figure B-3 The L band (21 cm) radiometer data (look angle = 0° and 40®) as 
a function of volumetric soil moisture content (%) in the 0-2.5 
layer 
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Figure B-4 The PMIS (2.8 cm) radiometer data (look angle = 49.2“) as a 

function of volumetric soil moisture content (%) in the 0-2.5 cm 
layer 




